Mouse RNA heatmap from human PC
Bing Shui
4/6/2020
Index
Library setup
suppressMessages(
suppressWarnings(
c(library(data.table),
library(gdata),
library(ggplot2),
library(RColorBrewer),
library(readr),
library(tidyverse),
library(plotly),
library(directlabels),
library(gplots),
library(knitr)
)))## [1] "data.table" "stats" "graphics" "grDevices" "utils"
## [6] "datasets" "methods" "base" "gdata" "data.table"
## [11] "stats" "graphics" "grDevices" "utils" "datasets"
## [16] "methods" "base" "ggplot2" "gdata" "data.table"
## [21] "stats" "graphics" "grDevices" "utils" "datasets"
## [26] "methods" "base" "RColorBrewer" "ggplot2" "gdata"
## [31] "data.table" "stats" "graphics" "grDevices" "utils"
## [36] "datasets" "methods" "base" "readr" "RColorBrewer"
## [41] "ggplot2" "gdata" "data.table" "stats" "graphics"
## [46] "grDevices" "utils" "datasets" "methods" "base"
## [51] "forcats" "stringr" "dplyr" "purrr" "tidyr"
## [56] "tibble" "tidyverse" "readr" "RColorBrewer" "ggplot2"
## [61] "gdata" "data.table" "stats" "graphics" "grDevices"
## [66] "utils" "datasets" "methods" "base" "plotly"
## [71] "forcats" "stringr" "dplyr" "purrr" "tidyr"
## [76] "tibble" "tidyverse" "readr" "RColorBrewer" "ggplot2"
## [81] "gdata" "data.table" "stats" "graphics" "grDevices"
## [86] "utils" "datasets" "methods" "base" "directlabels"
## [91] "plotly" "forcats" "stringr" "dplyr" "purrr"
## [96] "tidyr" "tibble" "tidyverse" "readr" "RColorBrewer"
## [101] "ggplot2" "gdata" "data.table" "stats" "graphics"
## [106] "grDevices" "utils" "datasets" "methods" "base"
## [111] "gplots" "directlabels" "plotly" "forcats" "stringr"
## [116] "dplyr" "purrr" "tidyr" "tibble" "tidyverse"
## [121] "readr" "RColorBrewer" "ggplot2" "gdata" "data.table"
## [126] "stats" "graphics" "grDevices" "utils" "datasets"
## [131] "methods" "base" "knitr" "gplots" "directlabels"
## [136] "plotly" "forcats" "stringr" "dplyr" "purrr"
## [141] "tidyr" "tibble" "tidyverse" "readr" "RColorBrewer"
## [146] "ggplot2" "gdata" "data.table" "stats" "graphics"
## [151] "grDevices" "utils" "datasets" "methods" "base"ID conversion from human gene symbol to mouse Ensembl ID
This conversion is done using ID database from MGI to be consistent with Doug Brubaker’s analysis using human RNA-Seq dataset.
Import the dataset containing Gene ID information from MGI
homolog_id <- read_csv("Human_2_mouse_symbol.csv")## Parsed with column specification:
## cols(
## `Human Marker Symbol` = col_character(),
## `Human Entrez Gene ID` = col_double(),
## `Mouse Marker Symbol` = col_character(),
## `MGI Marker Accession ID` = col_character()
## )ensembl_id <- read_csv("mouse_symbol_Ensembl.csv")## Parsed with column specification:
## cols(
## `MGI Accession ID` = col_character(),
## `Marker Symbol` = col_character(),
## `EntrezGene ID` = col_character(),
## `Ensembl Gene ID` = col_character(),
## `HomoloGene ID` = col_character()
## )# load the PC dataset
cCD_PC <- read_csv("cCD PC.csv")## Parsed with column specification:
## cols(
## Row = col_character(),
## cCD_PC1 = col_double(),
## cCD_PC2 = col_double(),
## cCD_PC3 = col_double(),
## cCD_PC4 = col_double(),
## cCD_PC5 = col_double(),
## cCD_PC6 = col_double(),
## cCD_PC7 = col_double(),
## cCD_PC8 = col_double(),
## cCD_PC9 = col_double(),
## cCD_PC10 = col_double(),
## cCD_PC11 = col_double(),
## cCD_PC12 = col_double(),
## cCD_PC13 = col_double(),
## cCD_PC14 = col_double(),
## cCD_PC15 = col_double(),
## cCD_PC16 = col_double()
## )iCD_PC <- read_csv("iCD PC.csv")## Parsed with column specification:
## cols(
## Row = col_character(),
## iCD_PC1 = col_double(),
## iCD_PC2 = col_double(),
## iCD_PC3 = col_double(),
## iCD_PC4 = col_double(),
## iCD_PC5 = col_double(),
## iCD_PC6 = col_double(),
## iCD_PC7 = col_double(),
## iCD_PC8 = col_double(),
## iCD_PC9 = col_double(),
## iCD_PC10 = col_double(),
## iCD_PC11 = col_double(),
## iCD_PC12 = col_double(),
## iCD_PC13 = col_double(),
## iCD_PC14 = col_double(),
## iCD_PC15 = col_double(),
## iCD_PC16 = col_double()
## )human_id <- as.data.frame(cCD_PC$Row)
human_id <- as_tibble(human_id)
colnames(human_id) <- "human_gene_symbol"
id_conversion <- inner_join(human_id, homolog_id, by = c("human_gene_symbol" = "Human Marker Symbol"))## Warning: Column `human_gene_symbol`/`Human Marker Symbol` joining factor and
## character vector, coercing into character vectornon_duplicates_idx <- which(duplicated(id_conversion$human_gene_symbol) == FALSE)
id_conversion <- id_conversion[non_duplicates_idx, ]
id_conversion <- inner_join(id_conversion, ensembl_id, by = c("Mouse Marker Symbol" = "Marker Symbol"))
non_duplicates_idx <- which(duplicated(id_conversion$`Mouse Marker Symbol`) == FALSE)
id_conversion <- id_conversion[non_duplicates_idx, ]
id_conversion <- id_conversion[,c(1,7)]
cCD_PC <- inner_join(cCD_PC, id_conversion, by = c("Row" = "human_gene_symbol"))
iCD_PC <- inner_join(iCD_PC, id_conversion, by = c("Row" = "human_gene_symbol"))Heatmap plotting based on cCD PC
For cCD PCs, PC1 is selected because it seperates diseases statuses. PC10 and PC14 are selected because they seperate TNF from TCT mice.
The goal of this is to examine how many of the genes that contribute to the seperation of TNF and TCT model based on human PCs are targeted by MK2 inhibitor.
For each PC, the top 325 contributing genes are selected for plotting the heatmap (10% of the list). Top contributing genes are genes with the largest absolute loading coefficients.
n = 10
threshold_p1 <- quantile(abs(cCD_PC$cCD_PC1), prob=1-n/100)
threshold_p10 <- quantile(abs(cCD_PC$cCD_PC10), prob=1-n/100)
threshold_p14 <- quantile(abs(cCD_PC$cCD_PC14), prob=1-n/100)
ensembl_p1 <- cCD_PC$`Ensembl Gene ID`[abs(cCD_PC$cCD_PC1) > threshold_p1]
ensembl_p10 <- cCD_PC$`Ensembl Gene ID`[abs(cCD_PC$cCD_PC10) > threshold_p10]
ensembl_p14 <- cCD_PC$`Ensembl Gene ID`[abs(cCD_PC$cCD_PC14) > threshold_p14]cCD_PC1
I will use I_V vs NI_V data to plot this heatmap since this is the PC that separates diseases status.
TCT_iv_niv_vst <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TCT/VST_I_V vs NI_V.csv")## Warning: Missing column names filled in: 'X1' [1]## Parsed with column specification:
## cols(
## X1 = col_character(),
## NI_V_1 = col_double(),
## NI_V_2 = col_double(),
## NI_V_4 = col_double(),
## NI_V_6 = col_double(),
## NI_V_7 = col_double(),
## I_V_1 = col_double(),
## I_V_2 = col_double(),
## I_V_3 = col_double(),
## I_V_4 = col_double(),
## I_V_5 = col_double(),
## I_V_6 = col_double()
## )TNF_iv_niv_vst <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TNF/VST_I_V vs NI_V.csv")## Warning: Missing column names filled in: 'X1' [1]## Parsed with column specification:
## cols(
## X1 = col_character(),
## NI_V_1 = col_double(),
## NI_V_2 = col_double(),
## NI_V_3 = col_double(),
## NI_V_4 = col_double(),
## NI_V_5 = col_double(),
## I_V_1 = col_double(),
## I_V_2 = col_double(),
## I_V_3 = col_double(),
## I_V_4 = col_double(),
## I_V_5 = col_double()
## )TCT_iv_niv_dif <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TCT/Differential Analysis_I_V vs NI_V.csv")## Warning: Missing column names filled in: 'X1' [1]## Parsed with column specification:
## cols(
## X1 = col_character(),
## baseMean = col_double(),
## log2FoldChange = col_double(),
## lfcSE = col_double(),
## stat = col_double(),
## pvalue = col_double(),
## padj = col_double()
## )TNF_iv_niv_dif <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TNF/Differential Analysis_I_V vs NI_V.csv")## Warning: Missing column names filled in: 'X1' [1]## Parsed with column specification:
## cols(
## X1 = col_character(),
## baseMean = col_double(),
## log2FoldChange = col_double(),
## lfcSE = col_double(),
## stat = col_double(),
## pvalue = col_double(),
## padj = col_double()
## )TCT model
All genes in the list of 325 from PC1
suppressMessages(library(mosaic))
pc1_genes <- TCT_iv_niv_vst[TCT_iv_niv_vst$X1 %in% ensembl_p1,]
pc1_genes <- as.data.frame(pc1_genes)
for (i in 1:dim(pc1_genes)[1]) {
pc1_genes[i,2:12] <- zscore(as.numeric(pc1_genes[i,2:12]))
}
my_palette <- colorRampPalette(c("blue", "white", "red"))(128)
heatmap_matrix <- as.matrix(pc1_genes[,2:12])
png('TCT_I_V vs NI_V in cCD PC1 Genes.png',
width = 600,
height = 1200,
res = 200,
pointsize = 8)
heatmap.2(heatmap_matrix,
main = "TCT I_V vs NI_V\nin cCD PC1 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2pdf('TCT_I_V vs NI_V in cCD PC1 Genes.pdf',
width = 6,
height = 10)
heatmap.2(heatmap_matrix,
main = "TCT I_V vs NI_V\nin cCD PC1 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2include_graphics("TCT_I_V vs NI_V in cCD PC1 Genes.png")
DEG from PC1
pc1_dif <- TCT_iv_niv_dif[TCT_iv_niv_dif$X1 %in% ensembl_p1,]
pc1_DEG <- subset(pc1_dif, pc1_dif$padj < 0.05 & abs(pc1_dif$log2FoldChange) > 0.25)
pc1_genes_DEG <- pc1_genes[pc1_genes$X1 %in% pc1_DEG$X1,]
pc1_genes_DEG <- inner_join(pc1_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
heatmap_matrix <- as.matrix(pc1_genes_DEG[,2:12])
png('TCT_I_V vs NI_V DEG in cCD PC1 Genes.png',
width = 800,
height = 600,
res = 200,
pointsize = 8)
heatmap.2(heatmap_matrix,
main = "TCT_I_V vs NI_V DEG\nin cCD PC1 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,3),
col=my_palette,
cexCol = 1,
margins = c(4,6),
trace = "none",
dendrogram = "row",
labRow = pc1_genes_DEG$`Marker Symbol`,
keysize = 2,
Colv = "NA")
dev.off()## quartz_off_screen
## 2pdf('TCT_I_V vs NI_V DEG in cCD PC1 Genes.pdf',
width = 8,
height = 6)
heatmap.2(heatmap_matrix,
main = "TCT_I_V vs NI_V DEG\nin cCD PC1 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,3),
col=my_palette,
cexCol = 1,
margins = c(4,6),
trace = "none",
dendrogram = "row",
labRow = pc1_genes_DEG$`Marker Symbol`,
keysize = 2,
Colv = "NA")
dev.off()## quartz_off_screen
## 2include_graphics("TCT_I_V vs NI_V DEG in cCD PC1 Genes.png")
6/324 genes are deemed significant.
TNF model
All genes in the list of 325 from PC1
pc1_genes <- TNF_iv_niv_vst[TNF_iv_niv_vst$X1 %in% ensembl_p1,]
pc1_genes <- as.data.frame(pc1_genes)
for (i in 1:dim(pc1_genes)[1]) {
pc1_genes[i,2:11] <- zscore(as.numeric(pc1_genes[i,2:11]))
}
heatmap_matrix <- as.matrix(pc1_genes[,2:11])
png('TNF_I_V vs NI_V in cCD PC1 Genes.png',
width = 600,
height = 1200,
res = 200,
pointsize = 8)
heatmap.2(heatmap_matrix,
main = "TNF I_V vs NI_V\nin cCD PC1 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2pdf('TNF_I_V vs NI_V in cCD PC1 Genes.pdf',
width = 6,
height = 10)
heatmap.2(heatmap_matrix,
main = "TNF I_V vs NI_V\nin cCD PC1 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2include_graphics("TNF_I_V vs NI_V in cCD PC1 Genes.png")
DEG from PC1
pc1_dif <- TNF_iv_niv_dif[TNF_iv_niv_dif$X1 %in% ensembl_p1,]
pc1_DEG <- subset(pc1_dif, pc1_dif$padj < 0.05 & abs(pc1_dif$log2FoldChange) > 0.25)
pc1_genes_DEG <- pc1_genes[pc1_genes$X1 %in% pc1_DEG$X1,]
pc1_genes_DEG <- inner_join(pc1_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
heatmap_matrix <- as.matrix(pc1_genes_DEG[,2:11])
png('TNF_I_V vs NI_V DEG in cCD PC1 Genes.png',
width = 600,
height = 1200,
res = 200,
pointsize = 8)
heatmap.2(heatmap_matrix,
main = "TNF\nI_V vs NI_V DEG\ncCD PC1 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(4,6),
trace = "none",
dendrogram = "row",
labRow = pc1_genes_DEG$`Marker Symbol`,
keysize = 2,
Colv = "NA")
dev.off()## quartz_off_screen
## 2pdf('TNF_I_V vs NI_V DEG in cCD PC1 Genes.pdf',
width = 6,
height = 10)
heatmap.2(heatmap_matrix,
main = "TNF\nI_V vs NI_V DEG\ncCD PC1 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(4,6),
trace = "none",
dendrogram = "row",
labRow = pc1_genes_DEG$`Marker Symbol`,
keysize = 2,
Colv = "NA")
dev.off()## quartz_off_screen
## 2include_graphics('TNF_I_V vs NI_V DEG in cCD PC1 Genes.png')
74/324 genes are deemed significant.
cCD_PC10 and cCD_PC14
I_MKI vs I_V datasets will be used here.
TCT_imki_iv_vst <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TCT/VST_I_MKI vs I_V.csv")## Warning: Missing column names filled in: 'X1' [1]## Parsed with column specification:
## cols(
## X1 = col_character(),
## I_V_1 = col_double(),
## I_V_2 = col_double(),
## I_V_3 = col_double(),
## I_V_4 = col_double(),
## I_V_5 = col_double(),
## I_V_6 = col_double(),
## I_MKI_1 = col_double(),
## I_MKI_2 = col_double(),
## I_MKI_3 = col_double(),
## I_MKI_4 = col_double(),
## I_MKI_5 = col_double(),
## I_MKI_6 = col_double()
## )TCT_nimki_niv_vst <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TCT/VST_NI_MKI vs NI_V.csv")## Warning: Missing column names filled in: 'X1' [1]## Parsed with column specification:
## cols(
## X1 = col_character(),
## NI_V_1 = col_double(),
## NI_V_2 = col_double(),
## NI_V_4 = col_double(),
## NI_V_6 = col_double(),
## NI_V_7 = col_double(),
## NI_MKI_1 = col_double(),
## NI_MKI_2 = col_double(),
## NI_MKI_3 = col_double(),
## NI_MKI_4 = col_double(),
## NI_MKI_5 = col_double(),
## NI_MKI_6 = col_double()
## )TCT_combined_vst <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TCT/VST_MKI vs V.csv")## Warning: Missing column names filled in: 'X1' [1]## Parsed with column specification:
## cols(
## .default = col_double(),
## X1 = col_character()
## )## See spec(...) for full column specifications.TNF_imki_iv_vst <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TNF/VST_I_MKI vs I_V.csv")## Warning: Missing column names filled in: 'X1' [1]## Parsed with column specification:
## cols(
## X1 = col_character(),
## I_V_1 = col_double(),
## I_V_2 = col_double(),
## I_V_3 = col_double(),
## I_V_4 = col_double(),
## I_V_5 = col_double(),
## I_MKI_1 = col_double(),
## I_MKI_2 = col_double(),
## I_MKI_3 = col_double(),
## I_MKI_4 = col_double(),
## I_MKI_5 = col_double(),
## I_MKI_6 = col_double()
## )TNF_nimki_niv_vst <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TNF/VST_NI_MKI vs NI_V.csv")## Warning: Missing column names filled in: 'X1' [1]## Parsed with column specification:
## cols(
## X1 = col_character(),
## NI_V_1 = col_double(),
## NI_V_2 = col_double(),
## NI_V_3 = col_double(),
## NI_V_4 = col_double(),
## NI_V_5 = col_double(),
## NI_MKI_1 = col_double(),
## NI_MKI_2 = col_double(),
## NI_MKI_3 = col_double(),
## NI_MKI_4 = col_double(),
## NI_MKI_5 = col_double()
## )TNF_combined_vst <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TNF/VST_MKI vs V.csv")## Warning: Missing column names filled in: 'X1' [1]## Parsed with column specification:
## cols(
## .default = col_double(),
## X1 = col_character()
## )
## See spec(...) for full column specifications.TCT_imki_iv_dif <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TCT/Differential Analysis_I_MKI vs I_V.csv")## Warning: Missing column names filled in: 'X1' [1]## Parsed with column specification:
## cols(
## X1 = col_character(),
## baseMean = col_double(),
## log2FoldChange = col_double(),
## lfcSE = col_double(),
## stat = col_double(),
## pvalue = col_double(),
## padj = col_double()
## )TCT_nimki_niv_dif <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TCT/Differential Analysis_NI_MKI vs NI_V.csv")## Warning: Missing column names filled in: 'X1' [1]## Parsed with column specification:
## cols(
## X1 = col_character(),
## baseMean = col_double(),
## log2FoldChange = col_double(),
## lfcSE = col_double(),
## stat = col_double(),
## pvalue = col_double(),
## padj = col_double()
## )TCT_combined_dif <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TCT/Differential Analysis_MKI vs V.csv")## Warning: Missing column names filled in: 'X1' [1]## Parsed with column specification:
## cols(
## X1 = col_character(),
## baseMean = col_double(),
## log2FoldChange = col_double(),
## lfcSE = col_double(),
## stat = col_double(),
## pvalue = col_double(),
## padj = col_double()
## )TNF_imki_iv_dif <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TNF/Differential Analysis_I_MKI vs I_V.csv")## Warning: Missing column names filled in: 'X1' [1]## Parsed with column specification:
## cols(
## X1 = col_character(),
## baseMean = col_double(),
## log2FoldChange = col_double(),
## lfcSE = col_double(),
## stat = col_double(),
## pvalue = col_double(),
## padj = col_double()
## )TNF_nimki_niv_dif <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TNF/Differential Analysis_NI_MKI vs NI_V.csv")## Warning: Missing column names filled in: 'X1' [1]## Parsed with column specification:
## cols(
## X1 = col_character(),
## baseMean = col_double(),
## log2FoldChange = col_double(),
## lfcSE = col_double(),
## stat = col_double(),
## pvalue = col_double(),
## padj = col_double()
## )TNF_combined_dif <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TNF/Differential Analysis_MKI vs V.csv")## Warning: Missing column names filled in: 'X1' [1]## Parsed with column specification:
## cols(
## X1 = col_character(),
## baseMean = col_double(),
## log2FoldChange = col_double(),
## lfcSE = col_double(),
## stat = col_double(),
## pvalue = col_double(),
## padj = col_double()
## )TCT model stratifying based on inflammation status
All genes in the list of 325 from PC10 and PC14
pc10_genes <- as.data.frame(TCT_imki_iv_vst[TCT_imki_iv_vst$X1 %in% ensembl_p10,])
pc14_genes <- as.data.frame(TCT_imki_iv_vst[TCT_imki_iv_vst$X1 %in% ensembl_p14,])
for (i in 1:dim(pc10_genes)[1]) {
pc10_genes[i,2:13] <- zscore(as.numeric(pc10_genes[i,2:13]))
}
for (i in 1:dim(pc14_genes)[1]) {
pc14_genes[i,2:13] <- zscore(as.numeric(pc14_genes[i,2:13]))
}
heatmap_matrix_pc10 <- as.matrix(pc10_genes[,2:13])
heatmap_matrix_pc14 <- as.matrix(pc14_genes[,2:13])
png('TCT_I_MKI vs I_V in cCD PC10 Genes.png',
width = 600,
height = 1200,
res = 200,
pointsize = 8)
heatmap.2(heatmap_matrix_pc10,
main = "TCT I_MKI vs I_V\nin cCD PC10 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2pdf('TCT_I_MKI vs I_V in cCD PC10 Genes.pdf',
width = 6,
height = 10)
heatmap.2(heatmap_matrix_pc10,
main = "TCT I_MKI vs I_V\nin cCD PC10 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2include_graphics('TCT_I_MKI vs I_V in cCD PC10 Genes.png')
png('TCT_I_MKI vs I_V in cCD PC14 Genes.png',
width = 600,
height = 1200,
res = 200,
pointsize = 8)
heatmap.2(heatmap_matrix_pc14,
main = "TCT I_MKI vs I_V\nin cCD PC14 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2pdf('TCT_I_MKI vs I_V in cCD PC14 Genes.pdf',
width = 6,
height = 10)
heatmap.2(heatmap_matrix_pc14,
main = "TCT I_MKI vs I_V\nin cCD PC14 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2include_graphics('TCT_I_MKI vs I_V in cCD PC14 Genes.png')
DEG from PC10 and PC14. The standard of DEG is relaxed to p < 0.1 here to keep consistent with previous analysis and LFC > 0.1.
pc10_dif <- TCT_imki_iv_dif[TCT_imki_iv_dif$X1 %in% ensembl_p10,]
pc14_dif <- TCT_imki_iv_dif[TCT_imki_iv_dif$X1 %in% ensembl_p14,]
pc10_DEG <- subset(pc10_dif, pc10_dif$padj < 0.1 & abs(pc10_dif$log2FoldChange > 0.1))
pc14_DEG <- subset(pc14_dif, pc14_dif$padj < 0.1 & abs(pc14_dif$log2FoldChange > 0.1))
pc10_genes_DEG <- pc10_genes[pc10_genes$X1 %in% pc10_DEG$X1,]
pc14_genes_DEG <- pc14_genes[pc14_genes$X1 %in% pc14_DEG$X1,]
pc10_genes_DEG <- inner_join(pc10_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
pc14_genes_DEG <- inner_join(pc14_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
print(paste("The number of MK2 inhibitor targets from TCT model from cCD PC10:", dim(pc10_genes_DEG)[1]))## [1] "The number of MK2 inhibitor targets from TCT model from cCD PC10: 0"print(paste("The number of MK2 inhibitor targets from TCT model from cCD PC14:", dim(pc14_genes_DEG)[1]))## [1] "The number of MK2 inhibitor targets from TCT model from cCD PC14: 1"We cannot make heatmap with only 1 item as input. So I will just print out the item here.
pc14_DEG## # A tibble: 1 x 7
## X1 baseMean log2FoldChange lfcSE stat pvalue padj
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ENSMUSG00000050022 55.3 1.34 0.263 4.82 0.00000141 0.000974I would also like to plot the Loading plot for this gene in PC10 and PC14 space.
loading_cord <- cCD_PC[cCD_PC$`Ensembl Gene ID` %in% pc14_genes_DEG$X1,]
loading_cord <- loading_cord[, c('cCD_PC10', 'cCD_PC14', 'Ensembl Gene ID')]
loading_cord <- inner_join(loading_cord, ensembl_id[,c(2,4)], by = c("Ensembl Gene ID"="Ensembl Gene ID"))
loading_cord$`Ensembl Gene ID` = NULL
loading_cord[2,] <- loading_cord[1,]
loading_cord <- as.data.frame(loading_cord)
loading_cord[1,] <- c(0,0,NA)
loading_cord$cCD_PC10 <- as.numeric(loading_cord$cCD_PC10)
loading_cord$cCD_PC14 <- as.numeric(loading_cord$cCD_PC14)
ggplot(loading_cord, aes(x=cCD_PC14, y=cCD_PC10)) +
geom_line(color = "red",
arrow = arrow()) +
xlim(-0.1, 0.1) +
ylim(-0.1, 0.1) +
ggtitle("Loading plot for TCT MK2i target in cCD PC10 and PC14 space") +
geom_dl(aes(label = loading_cord$`Marker Symbol`), method = list(dl.trans(x = x + 0.2), "last.points", cex = 0.8))## Warning: Use of `loading_cord$`Marker Symbol`` is discouraged. Use `Marker
## Symbol` instead.
TCT model disregarding inflammation status
All genes in the list of 325 from PC10 and PC14
pc10_genes <- as.data.frame(TCT_combined_vst[TCT_combined_vst$X1 %in% ensembl_p10,])
pc14_genes <- as.data.frame(TCT_combined_vst[TCT_combined_vst$X1 %in% ensembl_p14,])
for (i in 1:dim(pc10_genes)[1]) {
pc10_genes[i,2:26] <- zscore(as.numeric(pc10_genes[i,2:26]))
}
for (i in 1:dim(pc14_genes)[1]) {
pc14_genes[i,2:26] <- zscore(as.numeric(pc14_genes[i,2:26]))
}
heatmap_matrix_pc10 <- as.matrix(pc10_genes[,2:26])
heatmap_matrix_pc14 <- as.matrix(pc14_genes[,2:26])
png('TCT_MKI vs V in cCD PC10 Genes.png',
width = 800,
height = 1200,
res = 200,
pointsize = 8)
heatmap.2(heatmap_matrix_pc10,
main = "TCT MKI vs V\nin cCD PC10 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2pdf('TCT_MKI vs V in cCD PC10 Genes.pdf',
width = 8,
height = 10)
heatmap.2(heatmap_matrix_pc10,
main = "TCT MKI vs V\nin cCD PC10 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2include_graphics('TCT_MKI vs V in cCD PC10 Genes.png')
png('TCT_MKI vs V in cCD PC14 Genes.png',
width = 800,
height = 1200,
res = 200,
pointsize = 8)
heatmap.2(heatmap_matrix_pc14,
main = "TCT MKI vs V\nin cCD PC14 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2pdf('TCT_MKI vs V in cCD PC14 Genes.pdf',
width = 8,
height = 10)
heatmap.2(heatmap_matrix_pc14,
main = "TCT MKI vs V\nin cCD PC14 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2include_graphics('TCT_MKI vs V in cCD PC14 Genes.png')
DEG from PC10 and PC14. The standard of DEG is relaxed to p < 0.1 here to keep consistent with previous analysis and LFC > 0.1.
pc10_dif <- TCT_combined_dif[TCT_combined_dif$X1 %in% ensembl_p10,]
pc14_dif <- TCT_combined_dif[TCT_combined_dif$X1 %in% ensembl_p14,]
pc10_DEG <- subset(pc10_dif, pc10_dif$padj < 0.1 & abs(pc10_dif$log2FoldChange) > 0.1)
pc14_DEG <- subset(pc14_dif, pc14_dif$padj < 0.1 & abs(pc14_dif$log2FoldChange) > 0.1)
pc10_genes_DEG <- pc10_genes[pc10_genes$X1 %in% pc10_DEG$X1,]
pc14_genes_DEG <- pc14_genes[pc14_genes$X1 %in% pc14_DEG$X1,]
pc10_genes_DEG <- inner_join(pc10_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
pc14_genes_DEG <- inner_join(pc14_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
print(paste("The number of MK2 inhibitor targets from combined TCT model from cCD PC10:", dim(pc10_genes_DEG)[1]))## [1] "The number of MK2 inhibitor targets from combined TCT model from cCD PC10: 0"print(paste("The number of MK2 inhibitor targets from combined TCT model from cCD PC14:", dim(pc14_genes_DEG)[1]))## [1] "The number of MK2 inhibitor targets from combined TCT model from cCD PC14: 0"TNF model stratifying based on inflammation status
All genes in the list of 325 from PC10 and PC14
pc10_genes <- as.data.frame(TNF_imki_iv_vst[TNF_imki_iv_vst$X1 %in% ensembl_p10,])
pc14_genes <- as.data.frame(TNF_imki_iv_vst[TNF_imki_iv_vst$X1 %in% ensembl_p14,])
for (i in 1:dim(pc10_genes)[1]) {
pc10_genes[i,2:12] <- zscore(as.numeric(pc10_genes[i,2:12]))
}
for (i in 1:dim(pc14_genes)[1]) {
pc14_genes[i,2:12] <- zscore(as.numeric(pc14_genes[i,2:12]))
}
heatmap_matrix_pc10 <- as.matrix(pc10_genes[,2:12])
heatmap_matrix_pc14 <- as.matrix(pc14_genes[,2:12])
png('TNF_I_MKI vs I_V in cCD PC10 Genes.png',
width = 600,
height = 1200,
res = 200,
pointsize = 8)
heatmap.2(heatmap_matrix_pc10,
main = "TNF I_MKI vs I_V\nin cCD PC10 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2pdf('TNF_I_MKI vs I_V in cCD PC10 Genes.pdf',
width = 6,
height = 10)
heatmap.2(heatmap_matrix_pc10,
main = "TNF I_MKI vs I_V\nin cCD PC10 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2include_graphics('TNF_I_MKI vs I_V in cCD PC10 Genes.png')
png('TNF_I_MKI vs I_V in cCD PC14 Genes.png',
width = 600,
height = 1200,
res = 200,
pointsize = 8)
heatmap.2(heatmap_matrix_pc14,
main = "TNF I_MKI vs I_V\nin cCD PC14 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2pdf('TNF_I_MKI vs I_V in cCD PC14 Genes.pdf',
width = 6,
height = 10)
heatmap.2(heatmap_matrix_pc14,
main = "TNF I_MKI vs I_V\nin cCD PC14 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2include_graphics('TNF_I_MKI vs I_V in cCD PC14 Genes.png')
DEG from PC10 and PC14. The standard of DEG is relaxed to p < 0.1 here to keep consistent with previous analysis and LFC > 0.1.
pc10_dif <- TNF_imki_iv_dif[TNF_imki_iv_dif$X1 %in% ensembl_p10,]
pc14_dif <- TNF_imki_iv_dif[TNF_imki_iv_dif$X1 %in% ensembl_p14,]
pc10_DEG <- subset(pc10_dif, pc10_dif$padj < 0.1 & abs(pc10_dif$log2FoldChange > 0.1))
pc14_DEG <- subset(pc14_dif, pc14_dif$padj < 0.1 & abs(pc14_dif$log2FoldChange > 0.1))
pc10_genes_DEG <- pc10_genes[pc10_genes$X1 %in% pc10_DEG$X1,]
pc14_genes_DEG <- pc14_genes[pc14_genes$X1 %in% pc14_DEG$X1,]
pc10_genes_DEG <- inner_join(pc10_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
pc14_genes_DEG <- inner_join(pc14_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
print(paste("The number of MK2 inhibitor targets from TNF model from cCD PC10:", dim(pc10_genes_DEG)[1]))## [1] "The number of MK2 inhibitor targets from TNF model from cCD PC10: 0"print(paste("The number of MK2 inhibitor targets from TNF model from cCD PC14:", dim(pc14_genes_DEG)[1]))## [1] "The number of MK2 inhibitor targets from TNF model from cCD PC14: 0"TNF model disregarding inflammation status
All genes in the list of 325 from PC10 and PC14
pc10_genes <- as.data.frame(TNF_combined_vst[TNF_combined_vst$X1 %in% ensembl_p10,])
pc14_genes <- as.data.frame(TNF_combined_vst[TNF_combined_vst$X1 %in% ensembl_p14,])
for (i in 1:dim(pc10_genes)[1]) {
pc10_genes[i,2:22] <- zscore(as.numeric(pc10_genes[i,2:22]))
}
for (i in 1:dim(pc14_genes)[1]) {
pc14_genes[i,2:22] <- zscore(as.numeric(pc14_genes[i,2:22]))
}
heatmap_matrix_pc10 <- as.matrix(pc10_genes[,2:22])
heatmap_matrix_pc14 <- as.matrix(pc14_genes[,2:22])
png('TNF_MKI vs V in cCD PC10 Genes.png',
width = 800,
height = 1200,
res = 200,
pointsize = 8)
heatmap.2(heatmap_matrix_pc10,
main = "TNF MKI vs V\nin cCD PC10 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2pdf('TNF_MKI vs V in cCD PC10 Genes.pdf',
width = 8,
height = 10)
heatmap.2(heatmap_matrix_pc10,
main = "TNF MKI vs V\nin cCD PC10 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2include_graphics('TNF_MKI vs V in cCD PC10 Genes.png')
png('TNF_MKI vs V in cCD PC14 Genes.png',
width = 800,
height = 1200,
res = 200,
pointsize = 8)
heatmap.2(heatmap_matrix_pc14,
main = "TNF MKI vs V\nin cCD PC14 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2pdf('TNF_MKI vs V in cCD PC14 Genes.pdf',
width = 8,
height = 10)
heatmap.2(heatmap_matrix_pc14,
main = "TNF MKI vs V\nin cCD PC14 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2include_graphics('TNF_MKI vs V in cCD PC14 Genes.png')
DEG from PC10 and PC14. The standard of DEG is relaxed to p < 0.1 here to keep consistent with previous analysis and LFC > 0.1.
pc10_dif <- TNF_combined_dif[TNF_combined_dif$X1 %in% ensembl_p10,]
pc14_dif <- TNF_combined_dif[TNF_combined_dif$X1 %in% ensembl_p14,]
pc10_DEG <- subset(pc10_dif, pc10_dif$padj < 0.1 & abs(pc10_dif$log2FoldChange) > 0.1)
pc14_DEG <- subset(pc14_dif, pc14_dif$padj < 0.1 & abs(pc14_dif$log2FoldChange) > 0.1)
pc10_genes_DEG <- pc10_genes[pc10_genes$X1 %in% pc10_DEG$X1,]
pc14_genes_DEG <- pc14_genes[pc14_genes$X1 %in% pc14_DEG$X1,]
pc10_genes_DEG <- inner_join(pc10_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
pc14_genes_DEG <- inner_join(pc14_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
print(paste("The number of MK2 inhibitor targets from TNF model from cCD PC10:", dim(pc10_genes_DEG)[1]))## [1] "The number of MK2 inhibitor targets from TNF model from cCD PC10: 1"print(paste("The number of MK2 inhibitor targets from TNF model from cCD PC14:", dim(pc14_genes_DEG)[1]))## [1] "The number of MK2 inhibitor targets from TNF model from cCD PC14: 0"We cannot make heatmap with only 1 item as input. So I will just print out the item here.
pc10_DEG## # A tibble: 1 x 7
## X1 baseMean log2FoldChange lfcSE stat pvalue padj
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ENSMUSG00000028358 102. 0.582 0.175 4.94 0.000000781 0.000874I would also like to plot the Loading plot for this gene in PC10 and PC14 space.
loading_cord <- cCD_PC[cCD_PC$`Ensembl Gene ID` %in% pc10_genes_DEG$X1,]
loading_cord <- loading_cord[, c('cCD_PC10', 'cCD_PC14', 'Ensembl Gene ID')]
loading_cord <- inner_join(loading_cord, ensembl_id[,c(2,4)], by = c("Ensembl Gene ID"="Ensembl Gene ID"))
loading_cord$`Ensembl Gene ID` = NULL
loading_cord[2,] <- loading_cord[1,]
loading_cord <- as.data.frame(loading_cord)
loading_cord[1,] <- c(0,0,NA)
loading_cord$cCD_PC10 <- as.numeric(loading_cord$cCD_PC10)
loading_cord$cCD_PC14 <- as.numeric(loading_cord$cCD_PC14)
ggplot(loading_cord, aes(x=cCD_PC14, y=cCD_PC10)) +
geom_line(color = "red",
arrow = arrow(ends = "first")) +
xlim(-0.1, 0.1) +
ylim(-0.1, 0.1) +
ggtitle("Loading plot for combined TNF MK2i target in cCD PC10 and PC14 space") +
geom_dl(aes(label = loading_cord$`Marker Symbol`), method = list(dl.trans(x = x + 0.2), "last.points", cex = 0.8))## Warning: Use of `loading_cord$`Marker Symbol`` is discouraged. Use `Marker
## Symbol` instead.
Heatmap plotting based on iCD PC
For iCD PCs, PC1, PC14 and PC16 are selected because they seperate TNF from TCT mice.
The goal of this is to examine how many of the genes that contribute to the seperation of TNF and TCT model based on human PCs are targeted by MK2 inhibitor.
For each PC, the top 325 contributing genes are selected for plotting the heatmap (10% of the list). Top contributing genes are genes with the largest absolute loading coefficients.
threshold_p1 <- quantile(abs(iCD_PC$iCD_PC1), prob=1-n/100)
threshold_p14 <- quantile(abs(iCD_PC$iCD_PC14), prob=1-n/100)
threshold_p16 <- quantile(abs(iCD_PC$iCD_PC16), prob=1-n/100)
ensembl_p1 <- iCD_PC$`Ensembl Gene ID`[abs(iCD_PC$iCD_PC1) > threshold_p1]
ensembl_p14 <- iCD_PC$`Ensembl Gene ID`[abs(iCD_PC$iCD_PC14) > threshold_p14]
ensembl_p16 <- iCD_PC$`Ensembl Gene ID`[abs(iCD_PC$iCD_PC16) > threshold_p16]I_MKI vs I_V datasets will be used here for the analysis.
TCT model stratifying based on inflammation status
All genes in the list of 325 from PC1, PC14 and PC16
pc1_genes <- as.data.frame(TCT_imki_iv_vst[TCT_imki_iv_vst$X1 %in% ensembl_p1,])
pc14_genes <- as.data.frame(TCT_imki_iv_vst[TCT_imki_iv_vst$X1 %in% ensembl_p14,])
pc16_genes <-as.data.frame(TCT_imki_iv_vst[TCT_imki_iv_vst$X1 %in% ensembl_p16,])
for (i in 1:dim(pc1_genes)[1]) {
pc1_genes[i,2:13] <- zscore(as.numeric(pc1_genes[i,2:13]))
}
for (i in 1:dim(pc14_genes)[1]) {
pc14_genes[i,2:13] <- zscore(as.numeric(pc14_genes[i,2:13]))
}
for (i in 1:dim(pc16_genes)[1]) {
pc16_genes[i,2:13] <- zscore(as.numeric(pc16_genes[i,2:13]))
}
heatmap_matrix_pc1 <- as.matrix(pc1_genes[,2:13])
heatmap_matrix_pc14 <- as.matrix(pc14_genes[,2:13])
heatmap_matrix_pc16 <- as.matrix(pc16_genes[,2:13])
png('TCT_I_MKI vs I_V in iCD PC1 Genes.png',
width = 600,
height = 1200,
res = 200,
pointsize = 8)
heatmap.2(heatmap_matrix_pc1,
main = "TCT I_MKI vs I_V\nin iCD PC1 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2pdf('TCT_I_MKI vs I_V in iCD PC1 Genes.pdf',
width = 6,
height = 10)
heatmap.2(heatmap_matrix_pc1,
main = "TCT I_MKI vs I_V\nin iCD PC1 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2include_graphics('TCT_I_MKI vs I_V in iCD PC1 Genes.png')
png('TCT_I_MKI vs I_V in iCD PC14 Genes.png',
width = 600,
height = 1200,
res = 200,
pointsize = 8)
heatmap.2(heatmap_matrix_pc14,
main = "TCT I_MKI vs I_V\nin iCD PC14 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2pdf('TCT_I_MKI vs I_V in iCD PC14 Genes.pdf',
width = 6,
height = 10)
heatmap.2(heatmap_matrix_pc14,
main = "TCT I_MKI vs I_V\nin iCD PC14 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2include_graphics('TCT_I_MKI vs I_V in iCD PC14 Genes.png')
png('TCT_I_MKI vs I_V in iCD PC16 Genes.png',
width = 600,
height = 1200,
res = 200,
pointsize = 8)
heatmap.2(heatmap_matrix_pc16,
main = "TCT I_MKI vs I_V\nin iCD PC16 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2pdf('TCT_I_MKI vs I_V in iCD PC16 Genes.pdf',
width = 6,
height = 10)
heatmap.2(heatmap_matrix_pc16,
main = "TCT I_MKI vs I_V\nin iCD PC16 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2include_graphics('TCT_I_MKI vs I_V in iCD PC16 Genes.png')
DEG from PC1, PC14 and PC16. The standard of DEG is relaxed to p < 0.1 here to keep consistent with previous analysis and LFC > 0.1.
pc1_dif <- TCT_imki_iv_dif[TCT_imki_iv_dif$X1 %in% ensembl_p1,]
pc14_dif <- TCT_imki_iv_dif[TCT_imki_iv_dif$X1 %in% ensembl_p14,]
pc16_dif <- TCT_imki_iv_dif[TCT_imki_iv_dif$X1 %in% ensembl_p16,]
pc1_DEG <- subset(pc1_dif, pc1_dif$padj < 0.1 & abs(pc1_dif$log2FoldChange > 0.1))
pc14_DEG <- subset(pc14_dif, pc14_dif$padj < 0.1 & abs(pc14_dif$log2FoldChange > 0.1))
pc16_DEG <- subset(pc16_dif, pc16_dif$padj < 0.1 & abs(pc16_dif$log2FoldChange > 0.1))
pc1_genes_DEG <- pc1_genes[pc1_genes$X1 %in% pc1_DEG$X1,]
pc14_genes_DEG <- pc14_genes[pc14_genes$X1 %in% pc14_DEG$X1,]
pc16_genes_DEG <- pc16_genes[pc16_genes$X1 %in% pc16_DEG$X1,]
pc1_genes_DEG <- inner_join(pc1_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
pc14_genes_DEG <- inner_join(pc14_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
pc16_genes_DEG <- inner_join(pc16_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
print(paste("The number of MK2 inhibitor targets from TCT model from iCD PC1:", dim(pc1_genes_DEG)[1]))## [1] "The number of MK2 inhibitor targets from TCT model from iCD PC1: 0"print(paste("The number of MK2 inhibitor targets from TCT model from iCD PC14:", dim(pc14_genes_DEG)[1]))## [1] "The number of MK2 inhibitor targets from TCT model from iCD PC14: 0"print(paste("The number of MK2 inhibitor targets from TCT model from iCD PC16:", dim(pc16_genes_DEG)[1]))## [1] "The number of MK2 inhibitor targets from TCT model from iCD PC16: 0"TCT model disregarding inflammation status
All genes in the list of 325 from PC1, PC14 and PC16
pc1_genes <- as.data.frame(TCT_combined_vst[TCT_combined_vst$X1 %in% ensembl_p1,])
pc14_genes <- as.data.frame(TCT_combined_vst[TCT_combined_vst$X1 %in% ensembl_p14,])
pc16_genes <-as.data.frame(TCT_combined_vst[TCT_combined_vst$X1 %in% ensembl_p16,])
for (i in 1:dim(pc1_genes)[1]) {
pc1_genes[i,2:26] <- zscore(as.numeric(pc1_genes[i,2:26]))
}
for (i in 1:dim(pc14_genes)[1]) {
pc14_genes[i,2:26] <- zscore(as.numeric(pc14_genes[i,2:26]))
}
for (i in 1:dim(pc16_genes)[1]) {
pc16_genes[i,2:26] <- zscore(as.numeric(pc16_genes[i,2:26]))
}
heatmap_matrix_pc1 <- as.matrix(pc1_genes[,2:26])
heatmap_matrix_pc14 <- as.matrix(pc14_genes[,2:26])
heatmap_matrix_pc16 <- as.matrix(pc16_genes[,2:26])
png('TCT_MKI vs V in iCD PC1 Genes.png',
width = 800,
height = 1200,
res = 200,
pointsize = 8)
heatmap.2(heatmap_matrix_pc1,
main = "TCT MKI vs V\nin iCD PC1 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2pdf('TCT_MKI vs V in iCD PC1 Genes.pdf',
width = 8,
height = 10)
heatmap.2(heatmap_matrix_pc1,
main = "TCT MKI vs V\nin iCD PC1 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2include_graphics("TCT_MKI vs V in iCD PC1 Genes.png")
png('TCT_MKI vs V in iCD PC14 Genes.png',
width = 800,
height = 1200,
res = 200,
pointsize = 8)
heatmap.2(heatmap_matrix_pc14,
main = "TCT MKI vs V\nin iCD PC14 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2pdf('TCT_MKI vs V in iCD PC14 Genes.pdf',
width = 8,
height = 10)
heatmap.2(heatmap_matrix_pc14,
main = "TCT MKI vs V\nin iCD PC14 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2include_graphics("TCT_MKI vs V in iCD PC14 Genes.png")
png('TCT_MKI vs V in iCD PC16 Genes.png',
width = 800,
height = 1200,
res = 200,
pointsize = 8)
heatmap.2(heatmap_matrix_pc16,
main = "TCT MKI vs V\nin iCD PC16 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2pdf('TCT_MKI vs V in iCD PC16 Genes.pdf',
width = 8,
height = 10)
heatmap.2(heatmap_matrix_pc16,
main = "TCT MKI vs V\nin iCD PC16 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2include_graphics("TCT_MKI vs V in iCD PC16 Genes.png")
DEG from PC1, PC14 and PC16. The standard of DEG is relaxed to p < 0.1 here to keep consistent with previous analysis and LFC > 0.1.
pc1_dif <- TCT_combined_dif[TCT_combined_dif$X1 %in% ensembl_p1,]
pc14_dif <- TCT_combined_dif[TCT_combined_dif$X1 %in% ensembl_p14,]
pc16_dif <- TCT_combined_dif[TCT_combined_dif$X1 %in% ensembl_p16,]
pc1_DEG <- subset(pc1_dif, pc1_dif$padj < 0.1 & abs(pc1_dif$log2FoldChange > 0.1))
pc14_DEG <- subset(pc14_dif, pc14_dif$padj < 0.1 & abs(pc14_dif$log2FoldChange > 0.1))
pc16_DEG <- subset(pc16_dif, pc16_dif$padj < 0.1 & abs(pc16_dif$log2FoldChange > 0.1))
pc1_genes_DEG <- pc1_genes[pc1_genes$X1 %in% pc1_DEG$X1,]
pc14_genes_DEG <- pc14_genes[pc14_genes$X1 %in% pc14_DEG$X1,]
pc16_genes_DEG <- pc16_genes[pc16_genes$X1 %in% pc16_DEG$X1,]
pc1_genes_DEG <- inner_join(pc1_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
pc14_genes_DEG <- inner_join(pc14_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
pc16_genes_DEG <- inner_join(pc16_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
print(paste("The number of MK2 inhibitor targets from TCT model from iCD PC1:", dim(pc1_genes_DEG)[1]))## [1] "The number of MK2 inhibitor targets from TCT model from iCD PC1: 0"print(paste("The number of MK2 inhibitor targets from TCT model from iCD PC14:", dim(pc14_genes_DEG)[1]))## [1] "The number of MK2 inhibitor targets from TCT model from iCD PC14: 0"print(paste("The number of MK2 inhibitor targets from TCT model from iCD PC16:", dim(pc16_genes_DEG)[1]))## [1] "The number of MK2 inhibitor targets from TCT model from iCD PC16: 0"TNF model stratifying based on inflammation status
All genes in the list of 325 from PC1, PC14 and PC16
pc1_genes <- as.data.frame(TNF_imki_iv_vst[TNF_imki_iv_vst$X1 %in% ensembl_p1,])
pc14_genes <- as.data.frame(TNF_imki_iv_vst[TNF_imki_iv_vst$X1 %in% ensembl_p14,])
pc16_genes <- as.data.frame(TNF_imki_iv_vst[TNF_imki_iv_vst$X1 %in% ensembl_p16,])
for (i in 1:dim(pc1_genes)[1]) {
pc1_genes[i,2:12] <- zscore(as.numeric(pc1_genes[i,2:12]))
}
for (i in 1:dim(pc14_genes)[1]) {
pc14_genes[i,2:12] <- zscore(as.numeric(pc14_genes[i,2:12]))
}
for (i in 1:dim(pc16_genes)[1]) {
pc16_genes[i,2:12] <- zscore(as.numeric(pc16_genes[i,2:12]))
}
heatmap_matrix_pc1 <- as.matrix(pc1_genes[,2:12])
heatmap_matrix_pc14 <- as.matrix(pc14_genes[,2:12])
heatmap_matrix_pc16 <- as.matrix(pc16_genes[,2:12])
png('TNF_I_MKI vs I_V in iCD PC1 Genes.png',
width = 600,
height = 1200,
res = 200,
pointsize = 8)
heatmap.2(heatmap_matrix_pc1,
main = "TNF I_MKI vs I_V\nin iCD PC1 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2pdf('TNF_I_MKI vs I_V in iCD PC1 Genes.pdf',
width = 6,
height = 10)
heatmap.2(heatmap_matrix_pc1,
main = "TNF I_MKI vs I_V\nin iCD PC1 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2include_graphics('TNF_I_MKI vs I_V in iCD PC1 Genes.png')
png('TNF_I_MKI vs I_V in iCD PC14 Genes.png',
width = 600,
height = 1200,
res = 200,
pointsize = 8)
heatmap.2(heatmap_matrix_pc14,
main = "TNF I_MKI vs I_V\nin iCD PC14 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2pdf('TNF_I_MKI vs I_V in iCD PC14 Genes.pdf',
width = 6,
height = 10)
heatmap.2(heatmap_matrix_pc14,
main = "TNF I_MKI vs I_V\nin iCD PC14 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2include_graphics('TNF_I_MKI vs I_V in iCD PC14 Genes.png')
png('TNF_I_MKI vs I_V in iCD PC16 Genes.png',
width = 600,
height = 1200,
res = 200,
pointsize = 8)
heatmap.2(heatmap_matrix_pc16,
main = "TNF I_MKI vs I_V\nin iCD PC16 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2pdf('TNF_I_MKI vs I_V in iCD PC16 Genes.pdf',
width = 6,
height = 10)
heatmap.2(heatmap_matrix_pc16,
main = "TNF I_MKI vs I_V\nin iCD PC16 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2include_graphics('TNF_I_MKI vs I_V in iCD PC16 Genes.png')
DEG from PC1, PC14 and PC16. The standard of DEG is relaxed to p < 0.1 here to keep consistent with previous analysis and LFC > 0.1.
pc1_dif <- TNF_imki_iv_dif[TNF_imki_iv_dif$X1 %in% ensembl_p1,]
pc14_dif <- TNF_imki_iv_dif[TNF_imki_iv_dif$X1 %in% ensembl_p14,]
pc16_dif <- TNF_imki_iv_dif[TNF_imki_iv_dif$X1 %in% ensembl_p16,]
pc1_DEG <- subset(pc1_dif, pc1_dif$padj < 0.1 & abs(pc1_dif$log2FoldChange > 0.1))
pc14_DEG <- subset(pc14_dif, pc14_dif$padj < 0.1 & abs(pc14_dif$log2FoldChange > 0.1))
pc16_DEG <- subset(pc16_dif, pc16_dif$padj < 0.1 & abs(pc16_dif$log2FoldChange > 0.1))
pc1_genes_DEG <- pc1_genes[pc1_genes$X1 %in% pc1_DEG$X1,]
pc14_genes_DEG <- pc14_genes[pc14_genes$X1 %in% pc14_DEG$X1,]
pc16_genes_DEG <- pc16_genes[pc16_genes$X1 %in% pc16_DEG$X1,]
pc1_genes_DEG <- inner_join(pc1_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
pc14_genes_DEG <- inner_join(pc14_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
pc16_genes_DEG <- inner_join(pc16_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
print(paste("The number of MK2 inhibitor targets from TNF model from iCD PC1:", dim(pc1_genes_DEG)[1]))## [1] "The number of MK2 inhibitor targets from TNF model from iCD PC1: 1"print(paste("The number of MK2 inhibitor targets from TNF model from iCD PC14:", dim(pc14_genes_DEG)[1]))## [1] "The number of MK2 inhibitor targets from TNF model from iCD PC14: 0"print(paste("The number of MK2 inhibitor targets from TNF model from iCD PC16:", dim(pc16_genes_DEG)[1]))## [1] "The number of MK2 inhibitor targets from TNF model from iCD PC16: 0"Let’s take a look at the only MK2i target from TNF model in iCD PC1 space.
pc1_DEG## # A tibble: 1 x 7
## X1 baseMean log2FoldChange lfcSE stat pvalue padj
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ENSMUSG00000030064 347. 0.818 0.193 4.33 0.0000146 0.0113I would also like to plot the Loading plot for this gene in PC1/PC14 space and PC1/PC16 space.
loading_cord <- iCD_PC[iCD_PC$`Ensembl Gene ID` %in% pc1_genes_DEG$X1,]
loading_cord <- loading_cord[, c('iCD_PC1', 'iCD_PC14', 'iCD_PC16', 'Ensembl Gene ID')]
loading_cord <- inner_join(loading_cord, ensembl_id[,c(2,4)], by = c("Ensembl Gene ID"="Ensembl Gene ID"))
loading_cord$`Ensembl Gene ID` = NULL
loading_cord[2,] <- loading_cord[1,]
loading_cord <- as.data.frame(loading_cord)
loading_cord[1,] <- c(0,0,0,NA)
loading_cord$iCD_PC1 <- as.numeric(loading_cord$iCD_PC1)
loading_cord$iCD_PC14 <- as.numeric(loading_cord$iCD_PC14)
loading_cord$iCD_PC16 <- as.numeric(loading_cord$iCD_PC16)
ggplot(loading_cord, aes(x=iCD_PC1, y=iCD_PC14)) +
geom_line(color = "red",
arrow = arrow(ends = "first")) +
xlim(-0.1, 0.1) +
ylim(-0.1, 0.1) +
ggtitle("Loading plot for TNF MK2i target in iCD PC1 and PC14 space") +
geom_dl(aes(label = loading_cord$`Marker Symbol`), method = list(dl.trans(x = x - 1), dl.trans(y = y - 0.5), "last.points", cex = 0.8))## Warning: Use of `loading_cord$`Marker Symbol`` is discouraged. Use `Marker
## Symbol` instead.
ggplot(loading_cord, aes(x=iCD_PC1, y=iCD_PC16)) +
geom_line(color = "red",
arrow = arrow(ends = "first")) +
xlim(-0.1, 0.1) +
ylim(-0.1, 0.1) +
ggtitle("Loading plot for TNF MK2i target in iCD PC1 and PC16 space") +
geom_dl(aes(label = loading_cord$`Marker Symbol`), method = list(dl.trans(x = x - 1), dl.trans(y = y - 0.5), "last.points", cex = 0.8))## Warning: Use of `loading_cord$`Marker Symbol`` is discouraged. Use `Marker
## Symbol` instead.
TNF model disregarding inflammation status
All genes in the list of 325 from PC1, PC14 and PC16
pc1_genes <- as.data.frame(TNF_combined_vst[TNF_combined_vst$X1 %in% ensembl_p1,])
pc14_genes <- as.data.frame(TNF_combined_vst[TNF_combined_vst$X1 %in% ensembl_p14,])
pc16_genes <-as.data.frame(TNF_combined_vst[TNF_combined_vst$X1 %in% ensembl_p16,])
for (i in 1:dim(pc1_genes)[1]) {
pc1_genes[i,2:22] <- zscore(as.numeric(pc1_genes[i,2:22]))
}
for (i in 1:dim(pc14_genes)[1]) {
pc14_genes[i,2:22] <- zscore(as.numeric(pc14_genes[i,2:22]))
}
for (i in 1:dim(pc16_genes)[1]) {
pc16_genes[i,2:22] <- zscore(as.numeric(pc16_genes[i,2:22]))
}
heatmap_matrix_pc1 <- as.matrix(pc1_genes[,2:22])
heatmap_matrix_pc14 <- as.matrix(pc14_genes[,2:22])
heatmap_matrix_pc16 <- as.matrix(pc16_genes[,2:22])
png('TNF_MKI vs V in iCD PC1 Genes.png',
width = 800,
height = 1200,
res = 200,
pointsize = 8)
heatmap.2(heatmap_matrix_pc1,
main = "TNF MKI vs V\nin iCD PC1 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2pdf('TNF_MKI vs V in iCD PC1 Genes.pdf',
width = 8,
height = 10)
heatmap.2(heatmap_matrix_pc1,
main = "TNF MKI vs V\nin iCD PC1 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2include_graphics('TNF_MKI vs V in iCD PC1 Genes.png')
png('TNF_MKI vs V in iCD PC14 Genes.png',
width = 800,
height = 1200,
res = 200,
pointsize = 8)
heatmap.2(heatmap_matrix_pc14,
main = "TNF MKI vs V\nin iCD PC14 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2pdf('TNF_MKI vs V in iCD PC14 Genes.pdf',
width = 8,
height = 10)
heatmap.2(heatmap_matrix_pc14,
main = "TNF MKI vs V\nin iCD PC14 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2include_graphics('TNF_MKI vs V in iCD PC14 Genes.png')
png('TNF_MKI vs V in iCD PC16 Genes.png',
width = 800,
height = 1200,
res = 200,
pointsize = 8)
heatmap.2(heatmap_matrix_pc16,
main = "TNF MKI vs V\nin iCD PC16 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2pdf('TNF_MKI vs V in iCD PC16 Genes.pdf',
width = 8,
height = 10)
heatmap.2(heatmap_matrix_pc16,
main = "TNF MKI vs V\nin iCD PC16 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2include_graphics('TNF_MKI vs V in iCD PC16 Genes.png')
DEG from PC1, PC14 and PC16. The standard of DEG is relaxed to p < 0.1 here to keep consistent with previous analysis and LFC > 0.1.
pc1_dif <- TNF_combined_dif[TNF_combined_dif$X1 %in% ensembl_p1,]
pc14_dif <- TNF_combined_dif[TNF_combined_dif$X1 %in% ensembl_p14,]
pc16_dif <- TNF_combined_dif[TNF_combined_dif$X1 %in% ensembl_p16,]
pc1_DEG <- subset(pc1_dif, pc1_dif$padj < 0.1 & abs(pc1_dif$log2FoldChange > 0.1))
pc14_DEG <- subset(pc14_dif, pc14_dif$padj < 0.1 & abs(pc14_dif$log2FoldChange > 0.1))
pc16_DEG <- subset(pc16_dif, pc16_dif$padj < 0.1 & abs(pc16_dif$log2FoldChange > 0.1))
pc1_genes_DEG <- pc1_genes[pc1_genes$X1 %in% pc1_DEG$X1,]
pc14_genes_DEG <- pc14_genes[pc14_genes$X1 %in% pc14_DEG$X1,]
pc16_genes_DEG <- pc16_genes[pc16_genes$X1 %in% pc16_DEG$X1,]
pc1_genes_DEG <- inner_join(pc1_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
pc14_genes_DEG <- inner_join(pc14_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
pc16_genes_DEG <- inner_join(pc16_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
print(paste("The number of MK2 inhibitor targets from TNF model from iCD PC1:", dim(pc1_genes_DEG)[1]))## [1] "The number of MK2 inhibitor targets from TNF model from iCD PC1: 0"print(paste("The number of MK2 inhibitor targets from TNF model from iCD PC14:", dim(pc14_genes_DEG)[1]))## [1] "The number of MK2 inhibitor targets from TNF model from iCD PC14: 1"print(paste("The number of MK2 inhibitor targets from TNF model from iCD PC16:", dim(pc16_genes_DEG)[1]))## [1] "The number of MK2 inhibitor targets from TNF model from iCD PC16: 0"Let’s take a look at the only MK2i target from TNF model in iCD PC1 space.
pc14_DEG## # A tibble: 1 x 7
## X1 baseMean log2FoldChange lfcSE stat pvalue padj
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ENSMUSG00000089824 248. 0.892 0.201 4.27 0.0000200 0.0212I would also like to plot the Loading plot for this gene in PC1/PC14 space and PC1/PC16 space.
loading_cord <- iCD_PC[iCD_PC$`Ensembl Gene ID` %in% pc14_genes_DEG$X1,]
loading_cord <- loading_cord[, c('iCD_PC1', 'iCD_PC14', 'iCD_PC16', 'Ensembl Gene ID')]
loading_cord <- inner_join(loading_cord, ensembl_id[,c(2,4)], by = c("Ensembl Gene ID"="Ensembl Gene ID"))
loading_cord$`Ensembl Gene ID` = NULL
loading_cord[2,] <- loading_cord[1,]
loading_cord <- as.data.frame(loading_cord)
loading_cord[1,] <- c(0,0,0,NA)
loading_cord$iCD_PC1 <- as.numeric(loading_cord$iCD_PC1)
loading_cord$iCD_PC14 <- as.numeric(loading_cord$iCD_PC14)
loading_cord$iCD_PC16 <- as.numeric(loading_cord$iCD_PC16)
ggplot(loading_cord, aes(x=iCD_PC1, y=iCD_PC14)) +
geom_line(color = "red",
arrow = arrow(ends = "first")) +
xlim(-0.1, 0.1) +
ylim(-0.1, 0.1) +
ggtitle("Loading plot for TNF MK2i target in iCD PC1 and PC14 space") +
geom_dl(aes(label = loading_cord$`Marker Symbol`), method = list(dl.trans(x = x - 1), dl.trans(y = y - 0.5), "last.points", cex = 0.8))## Warning: Use of `loading_cord$`Marker Symbol`` is discouraged. Use `Marker
## Symbol` instead.
ggplot(loading_cord, aes(x=iCD_PC1, y=iCD_PC16)) +
geom_line(color = "red",
arrow = arrow(ends = "first")) +
xlim(-0.1, 0.1) +
ylim(-0.1, 0.1) +
ggtitle("Loading plot for TNF MK2i target in iCD PC1 and PC16 space") +
geom_dl(aes(label = loading_cord$`Marker Symbol`), method = list(dl.trans(x = x - 1), dl.trans(y = y - 0.5), "last.points", cex = 0.8))## Warning: Use of `loading_cord$`Marker Symbol`` is discouraged. Use `Marker
## Symbol` instead.
Heatmap plotting based on UC PC using MK2i dataset
TCT and TNF MK2i treated dataset with and without inflammation are used here.
# ID conversion
# load the PC dataset
UC_PC <- read_csv("UC PC.csv")## Parsed with column specification:
## cols(
## Row = col_character(),
## Var1_1 = col_double(),
## Var1_2 = col_double(),
## Var1_3 = col_double(),
## Var1_4 = col_double(),
## Var1_5 = col_double(),
## Var1_6 = col_double(),
## Var1_7 = col_double(),
## Var1_8 = col_double(),
## Var1_9 = col_double(),
## Var1_10 = col_double()
## )human_id <- as.data.frame(UC_PC$Row)
human_id <- as_tibble(human_id)
colnames(human_id) <- "human_gene_symbol"
id_conversion <- inner_join(human_id, homolog_id, by = c("human_gene_symbol" = "Human Marker Symbol"))## Warning: Column `human_gene_symbol`/`Human Marker Symbol` joining factor and
## character vector, coercing into character vectornon_duplicates_idx <- which(duplicated(id_conversion$human_gene_symbol) == FALSE)
id_conversion <- id_conversion[non_duplicates_idx, ]
id_conversion <- inner_join(id_conversion, ensembl_id, by = c("Mouse Marker Symbol" = "Marker Symbol"))
non_duplicates_idx <- which(duplicated(id_conversion$`Mouse Marker Symbol`) == FALSE)
id_conversion <- id_conversion[non_duplicates_idx, ]
id_conversion <- id_conversion[,c(1,7)]
UC_PC <- inner_join(UC_PC, id_conversion, by = c("Row" = "human_gene_symbol"))For UC PCs, PC1 and PC2 are selected because they seperates MK2i treatment statuses in TCT mice.
The goal of this is to examine how many of the genes that contribute to the seperation of MK2i vs Vehicle treated TCT mice based on human PCs are targeted by MK2 inhibitor.
For each PC, the top 49 contributing genes are selected for plotting the heatmap (10% of the list). Top contributing genes are genes with the largest absolute loading coefficients.
n = 10
threshold_p1 <- quantile(abs(UC_PC$Var1_1), prob=1-n/100)
threshold_p2 <- quantile(abs(UC_PC$Var1_2), prob=1-n/100)
ensembl_p1 <- UC_PC$`Ensembl Gene ID`[abs(UC_PC$Var1_1) > threshold_p1]
ensembl_p2 <- UC_PC$`Ensembl Gene ID`[abs(UC_PC$Var1_2) > threshold_p2]TCT model with inflammation treated with MK2i
All genes in the list of 49 from PC1 and PC2
pc1_genes <- as.data.frame(TCT_imki_iv_vst[TCT_imki_iv_vst$X1 %in% ensembl_p1,])
pc2_genes <- as.data.frame(TCT_imki_iv_vst[TCT_imki_iv_vst$X1 %in% ensembl_p2,])
for (i in 1:dim(pc1_genes)[1]) {
pc1_genes[i,2:13] <- zscore(as.numeric(pc1_genes[i,2:13]))
}
for (i in 1:dim(pc2_genes)[1]) {
pc2_genes[i,2:13] <- zscore(as.numeric(pc2_genes[i,2:13]))
}
heatmap_matrix_pc1 <- as.matrix(pc1_genes[,2:13])
heatmap_matrix_pc2 <- as.matrix(pc2_genes[,2:13])
png('TCT_I_MKI vs I_V in UC PC1 Genes.png',
width = 600,
height = 1200,
res = 200,
pointsize = 8)
heatmap.2(heatmap_matrix_pc1,
main = "TCT I_MKI vs I_V\nin UC PC1 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2pdf('TCT_I_MKI vs I_V in UC PC1 Genes.pdf',
width = 6,
height = 10)
heatmap.2(heatmap_matrix_pc1,
main = "TCT I_MKI vs I_V\nin UC PC1 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2include_graphics('TCT_I_MKI vs I_V in UC PC1 Genes.png')
png('TCT_I_MKI vs I_V in UC PC2 Genes.png',
width = 600,
height = 1200,
res = 200,
pointsize = 8)
heatmap.2(heatmap_matrix_pc2,
main = "TCT I_MKI vs I_V\nin UC PC2 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2pdf('TCT_I_MKI vs I_V in UC PC2 Genes.pdf',
width = 6,
height = 10)
heatmap.2(heatmap_matrix_pc2,
main = "TCT I_MKI vs I_V\nin UC PC2 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2include_graphics('TCT_I_MKI vs I_V in UC PC2 Genes.png')
DEG from PC1 and PC2. The standard of DEG is relaxed to p < 0.1 here to keep consistent with previous analysis and LFC > 0.1.
pc1_dif <- TCT_imki_iv_dif[TCT_imki_iv_dif$X1 %in% ensembl_p1,]
pc2_dif <- TCT_imki_iv_dif[TCT_imki_iv_dif$X1 %in% ensembl_p2,]
pc1_dif## # A tibble: 49 x 7
## X1 baseMean log2FoldChange lfcSE stat pvalue padj
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ENSMUSG00000004994 84.8 -0.508 0.295 -1.57 0.115 1.00
## 2 ENSMUSG00000005299 1920. -0.119 0.184 -0.587 0.557 1.00
## 3 ENSMUSG00000006205 486. 0.318 0.257 1.09 0.274 1.00
## 4 ENSMUSG00000020212 64.7 -0.344 0.268 -1.22 0.223 1.00
## 5 ENSMUSG00000020541 839. -0.129 0.188 -0.613 0.540 1.00
## 6 ENSMUSG00000020719 8112. -0.220 0.247 -0.797 0.425 1.00
## 7 ENSMUSG00000020810 291. 0.194 0.195 0.883 0.377 1.00
## 8 ENSMUSG00000020865 1996. 0.482 0.263 1.52 0.128 1.00
## 9 ENSMUSG00000021356 44.2 0.343 0.290 0.879 0.379 1.00
## 10 ENSMUSG00000021428 224. -0.313 0.233 -1.19 0.234 1.00
## # … with 39 more rowspc2_dif## # A tibble: 49 x 7
## X1 baseMean log2FoldChange lfcSE stat pvalue padj
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ENSMUSG00000000441 1301. -0.364 0.147 -2.44 0.0148 1.00
## 2 ENSMUSG00000001870 236. 0.0587 0.232 0.0630 0.950 1.00
## 3 ENSMUSG00000008575 628. 0.262 0.236 0.918 0.359 1.00
## 4 ENSMUSG00000009555 1008. -0.139 0.169 -0.777 0.437 1.00
## 5 ENSMUSG00000010392 766. -0.308 0.244 -1.09 0.275 1.00
## 6 ENSMUSG00000012350 3401. -0.172 0.259 -0.528 0.597 1.00
## 7 ENSMUSG00000017418 579. -0.481 0.243 -1.87 0.0617 1.00
## 8 ENSMUSG00000017478 871. -0.351 0.278 -1.09 0.275 1.00
## 9 ENSMUSG00000019960 1159. -0.336 0.240 -1.30 0.192 1.00
## 10 ENSMUSG00000020149 5883. -0.163 0.131 -1.23 0.220 1.00
## # … with 39 more rowspc1_DEG <- subset(pc1_dif, pc1_dif$padj < 0.1 & abs(pc1_dif$log2FoldChange > 0.1))
pc2_DEG <- subset(pc2_dif, pc2_dif$padj < 0.1 & abs(pc2_dif$log2FoldChange > 0.1))
pc1_genes_DEG <- pc1_genes[pc1_genes$X1 %in% pc1_DEG$X1,]
pc2_genes_DEG <- pc2_genes[pc2_genes$X1 %in% pc2_DEG$X1,]
pc1_genes_DEG <- inner_join(pc1_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
pc2_genes_DEG <- inner_join(pc2_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
print(paste("The number of MK2 inhibitor targets from TCT model from UC PC1:", dim(pc1_genes_DEG)[1]))## [1] "The number of MK2 inhibitor targets from TCT model from UC PC1: 0"print(paste("The number of MK2 inhibitor targets from TCT model from UC PC2:", dim(pc2_genes_DEG)[1]))## [1] "The number of MK2 inhibitor targets from TCT model from UC PC2: 0"TCT model with no inflammation treated with MK2i
All genes in the list of 49 from PC1 and PC2
pc1_genes <- as.data.frame(TCT_nimki_niv_vst[TCT_nimki_niv_vst$X1 %in% ensembl_p1,])
pc2_genes <- as.data.frame(TCT_nimki_niv_vst[TCT_nimki_niv_vst$X1 %in% ensembl_p2,])
for (i in 1:dim(pc1_genes)[1]) {
pc1_genes[i,2:12] <- zscore(as.numeric(pc1_genes[i,2:12]))
}
for (i in 1:dim(pc2_genes)[1]) {
pc2_genes[i,2:12] <- zscore(as.numeric(pc2_genes[i,2:12]))
}
heatmap_matrix_pc1 <- as.matrix(pc1_genes[,2:12])
heatmap_matrix_pc2 <- as.matrix(pc2_genes[,2:12])
png('TCT_NI_MKI vs NI_V in UC PC1 Genes.png',
width = 600,
height = 1200,
res = 200,
pointsize = 8)
heatmap.2(heatmap_matrix_pc1,
main = "TCT NI_MKI vs NI_V\nin UC PC1 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2pdf('TCT_NI_MKI vs NI_V in UC PC1 Genes.pdf',
width = 6,
height = 10)
heatmap.2(heatmap_matrix_pc1,
main = "TCT NI_MKI vs NI_V\nin UC PC1 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2include_graphics('TCT_NI_MKI vs NI_V in UC PC1 Genes.png')
png('TCT_NI_MKI vs NI_V in UC PC2 Genes.png',
width = 600,
height = 1200,
res = 200,
pointsize = 8)
heatmap.2(heatmap_matrix_pc2,
main = "TCT NI_MKI vs NI_V\nin UC PC2 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2pdf('TCT_NI_MKI vs NI_V in UC PC2 Genes.pdf',
width = 6,
height = 10)
heatmap.2(heatmap_matrix_pc2,
main = "TCT NI_MKI vs NI_V\nin UC PC2 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2include_graphics('TCT_NI_MKI vs NI_V in UC PC2 Genes.png')
DEG from PC1 and PC2. The standard of DEG is relaxed to p < 0.1 here to keep consistent with previous analysis and LFC > 0.1.
pc1_dif <- TCT_nimki_niv_dif[TCT_nimki_niv_dif$X1 %in% ensembl_p1,]
pc2_dif <- TCT_nimki_niv_dif[TCT_nimki_niv_dif$X1 %in% ensembl_p2,]
pc1_dif## # A tibble: 49 x 7
## X1 baseMean log2FoldChange lfcSE stat pvalue padj
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ENSMUSG00000004994 102. -0.296 0.191 -1.65 0.0981 1.00
## 2 ENSMUSG00000005299 2443. -0.107 0.126 -0.823 0.410 1.00
## 3 ENSMUSG00000006205 507. 0.0214 0.188 0.0193 0.985 1.00
## 4 ENSMUSG00000020212 54.0 0.0328 0.194 0.118 0.906 1.00
## 5 ENSMUSG00000020541 1084. 0.131 0.101 1.32 0.186 1.00
## 6 ENSMUSG00000020719 9771. -0.0156 0.146 -0.0439 0.965 1.00
## 7 ENSMUSG00000020810 255. -0.144 0.166 -0.877 0.381 1.00
## 8 ENSMUSG00000020865 3098. -0.126 0.144 -0.892 0.372 1.00
## 9 ENSMUSG00000021356 12.2 -0.0571 0.172 -0.559 0.576 1.00
## 10 ENSMUSG00000021428 200. 0.0337 0.156 0.243 0.808 1.00
## # … with 39 more rowspc2_dif## # A tibble: 49 x 7
## X1 baseMean log2FoldChange lfcSE stat pvalue padj
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ENSMUSG00000000441 1486. 0.0277 0.142 0.159 0.873 1.00
## 2 ENSMUSG00000001870 319. 0.0423 0.161 0.183 0.855 1.00
## 3 ENSMUSG00000008575 1119. -0.0938 0.115 -0.768 0.443 1.00
## 4 ENSMUSG00000009555 999. 0.0492 0.142 0.334 0.739 1.00
## 5 ENSMUSG00000010392 1012. -0.0414 0.191 -0.0328 0.974 1.00
## 6 ENSMUSG00000012350 6681. -0.00547 0.165 0.0936 0.925 1.00
## 7 ENSMUSG00000017418 706. -0.135 0.175 -0.613 0.540 1.00
## 8 ENSMUSG00000017478 1033. 0.133 0.176 0.851 0.395 1.00
## 9 ENSMUSG00000019960 776. -0.00514 0.193 -0.0380 0.970 1.00
## 10 ENSMUSG00000020149 6450. -0.120 0.0941 -1.26 0.207 1.00
## # … with 39 more rowspc1_DEG <- subset(pc1_dif, pc1_dif$padj < 0.1 & abs(pc1_dif$log2FoldChange > 0.1))
pc2_DEG <- subset(pc2_dif, pc2_dif$padj < 0.1 & abs(pc2_dif$log2FoldChange > 0.1))
pc1_genes_DEG <- pc1_genes[pc1_genes$X1 %in% pc1_DEG$X1,]
pc2_genes_DEG <- pc2_genes[pc2_genes$X1 %in% pc2_DEG$X1,]
pc1_genes_DEG <- inner_join(pc1_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
pc2_genes_DEG <- inner_join(pc2_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
print(paste("The number of MK2 inhibitor targets from TCT model from UC PC1:", dim(pc1_genes_DEG)[1]))## [1] "The number of MK2 inhibitor targets from TCT model from UC PC1: 0"print(paste("The number of MK2 inhibitor targets from TCT model from UC PC2:", dim(pc2_genes_DEG)[1]))## [1] "The number of MK2 inhibitor targets from TCT model from UC PC2: 0"TNF model with inflammation treated with MK2i
All genes in the list of 49 from PC1 and PC2
pc1_genes <- as.data.frame(TNF_imki_iv_vst[TNF_imki_iv_vst$X1 %in% ensembl_p1,])
pc2_genes <- as.data.frame(TNF_imki_iv_vst[TNF_imki_iv_vst$X1 %in% ensembl_p2,])
for (i in 1:dim(pc1_genes)[1]) {
pc1_genes[i,2:12] <- zscore(as.numeric(pc1_genes[i,2:12]))
}
for (i in 1:dim(pc2_genes)[1]) {
pc2_genes[i,2:12] <- zscore(as.numeric(pc2_genes[i,2:12]))
}
heatmap_matrix_pc1 <- as.matrix(pc1_genes[,2:12])
heatmap_matrix_pc2 <- as.matrix(pc2_genes[,2:12])
png('TNF_I_MKI vs I_V in UC PC1 Genes.png',
width = 600,
height = 1200,
res = 200,
pointsize = 8)
heatmap.2(heatmap_matrix_pc1,
main = "TNF I_MKI vs I_V\nin UC PC1 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2pdf('TNF_I_MKI vs I_V in UC PC1 Genes.pdf',
width = 6,
height = 10)
heatmap.2(heatmap_matrix_pc1,
main = "TNF I_MKI vs I_V\nin UC PC1 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2include_graphics('TNF_I_MKI vs I_V in UC PC1 Genes.png')
png('TNF_I_MKI vs I_V in UC PC2 Genes.png',
width = 600,
height = 1200,
res = 200,
pointsize = 8)
heatmap.2(heatmap_matrix_pc2,
main = "TNF I_MKI vs I_V\nin UC PC2 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2pdf('TNF_I_MKI vs I_V in UC PC2 Genes.pdf',
width = 6,
height = 10)
heatmap.2(heatmap_matrix_pc2,
main = "TNF I_MKI vs I_V\nin UC PC2 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2include_graphics('TNF_I_MKI vs I_V in UC PC2 Genes.png')
DEG from PC1 and PC2. The standard of DEG is relaxed to p < 0.1 here to keep consistent with previous analysis and LFC > 0.1.
pc1_dif <- TNF_imki_iv_dif[TNF_imki_iv_dif$X1 %in% ensembl_p1,]
pc2_dif <- TNF_imki_iv_dif[TNF_imki_iv_dif$X1 %in% ensembl_p2,]
pc1_dif## # A tibble: 49 x 7
## X1 baseMean log2FoldChange lfcSE stat pvalue padj
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ENSMUSG00000004994 48.9 0.150 0.246 0.603 0.547 1.00
## 2 ENSMUSG00000005299 1109. -0.126 0.108 -1.14 0.255 1.00
## 3 ENSMUSG00000006205 470. 0.0540 0.240 -0.148 0.883 1.00
## 4 ENSMUSG00000020212 31.5 -0.0306 0.239 -0.217 0.828 1.00
## 5 ENSMUSG00000020541 503. -0.0682 0.157 -0.496 0.620 1.00
## 6 ENSMUSG00000020719 5321. -0.216 0.119 -1.79 0.0738 1.00
## 7 ENSMUSG00000020810 148. 0.366 0.212 1.62 0.105 1.00
## 8 ENSMUSG00000020865 903. -0.0821 0.168 -0.339 0.735 1.00
## 9 ENSMUSG00000021356 242. 0.162 0.203 0.558 0.577 1.00
## 10 ENSMUSG00000021428 119. -0.111 0.187 -0.455 0.649 1.00
## # … with 39 more rowspc2_dif## # A tibble: 49 x 7
## X1 baseMean log2FoldChange lfcSE stat pvalue padj
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ENSMUSG00000000441 751. -0.164 0.142 -1.08 0.278 1.00
## 2 ENSMUSG00000001870 232. 0.375 0.235 1.22 0.224 1.00
## 3 ENSMUSG00000008575 442. 0.0164 0.202 0.0658 0.948 1.00
## 4 ENSMUSG00000009555 575. 0.264 0.184 1.47 0.140 1.00
## 5 ENSMUSG00000010392 499. 0.0190 0.151 -0.122 0.903 1.00
## 6 ENSMUSG00000012350 967. 0.382 0.209 1.66 0.0960 1.00
## 7 ENSMUSG00000017418 311. 0.0123 0.207 0.0503 0.960 1.00
## 8 ENSMUSG00000017478 538. 0.210 0.242 0.949 0.342 1.00
## 9 ENSMUSG00000019960 584. -0.287 0.193 -1.32 0.187 1.00
## 10 ENSMUSG00000020149 3198. -0.127 0.153 -0.873 0.383 1.00
## # … with 39 more rowspc1_DEG <- subset(pc1_dif, pc1_dif$padj < 0.1 & abs(pc1_dif$log2FoldChange > 0.1))
pc2_DEG <- subset(pc2_dif, pc2_dif$padj < 0.1 & abs(pc2_dif$log2FoldChange > 0.1))
pc1_genes_DEG <- pc1_genes[pc1_genes$X1 %in% pc1_DEG$X1,]
pc2_genes_DEG <- pc2_genes[pc2_genes$X1 %in% pc2_DEG$X1,]
pc1_genes_DEG <- inner_join(pc1_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
pc2_genes_DEG <- inner_join(pc2_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
print(paste("The number of MK2 inhibitor targets from TNF model from UC PC1:", dim(pc1_genes_DEG)[1]))## [1] "The number of MK2 inhibitor targets from TNF model from UC PC1: 0"print(paste("The number of MK2 inhibitor targets from TNF model from UC PC2:", dim(pc2_genes_DEG)[1]))## [1] "The number of MK2 inhibitor targets from TNF model from UC PC2: 0"TNF model with no inflammation treated with MK2i
All genes in the list of 49 from PC1 and PC2
pc1_genes <- as.data.frame(TNF_nimki_niv_vst[TNF_nimki_niv_vst$X1 %in% ensembl_p1,])
pc2_genes <- as.data.frame(TNF_nimki_niv_vst[TNF_nimki_niv_vst$X1 %in% ensembl_p2,])
for (i in 1:dim(pc1_genes)[1]) {
pc1_genes[i,2:11] <- zscore(as.numeric(pc1_genes[i,2:11]))
}
for (i in 1:dim(pc2_genes)[1]) {
pc2_genes[i,2:11] <- zscore(as.numeric(pc2_genes[i,2:11]))
}
heatmap_matrix_pc1 <- as.matrix(pc1_genes[,2:11])
heatmap_matrix_pc2 <- as.matrix(pc2_genes[,2:11])
png('TNF_NI_MKI vs NI_V in UC PC1 Genes.png',
width = 600,
height = 1200,
res = 200,
pointsize = 8)
heatmap.2(heatmap_matrix_pc1,
main = "TNF NI_MKI vs NI_V\nin UC PC1 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2pdf('TNF_NI_MKI vs NI_V in UC PC1 Genes.pdf',
width = 6,
height = 10)
heatmap.2(heatmap_matrix_pc1,
main = "TNF NI_MKI vs NI_V\nin UC PC1 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2include_graphics('TNF_NI_MKI vs NI_V in UC PC1 Genes.png')
png('TNF_NI_MKI vs NI_V in UC PC2 Genes.png',
width = 600,
height = 1200,
res = 200,
pointsize = 8)
heatmap.2(heatmap_matrix_pc2,
main = "TNF NI_MKI vs NI_V\nin UC PC2 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2pdf('TNF_NI_MKI vs NI_V in UC PC2 Genes.pdf',
width = 6,
height = 10)
heatmap.2(heatmap_matrix_pc2,
main = "TNF NI_MKI vs NI_V\nin UC PC2 Genes",
density.info = "none",
key = TRUE,
lhei = c(1,7),
col=my_palette,
cexCol = 1,
margins = c(8,2),
trace = "none",
dendrogram = "row",
labRow = FALSE,
keysize = 2,
ylab = "Genes",
Colv = "NA")
dev.off()## quartz_off_screen
## 2include_graphics('TNF_NI_MKI vs NI_V in UC PC2 Genes.png')
DEG from PC1 and PC2. The standard of DEG is relaxed to p < 0.1 here to keep consistent with previous analysis and LFC > 0.1.
pc1_dif <- TNF_nimki_niv_dif[TNF_nimki_niv_dif$X1 %in% ensembl_p1,]
pc2_dif <- TNF_nimki_niv_dif[TNF_nimki_niv_dif$X1 %in% ensembl_p2,]
pc1_dif## # A tibble: 49 x 7
## X1 baseMean log2FoldChange lfcSE stat pvalue padj
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ENSMUSG00000004994 73.1 0.0694 0.198 0.216 0.829 1.00
## 2 ENSMUSG00000005299 1924. 0.00207 0.146 0.0132 0.989 1.00
## 3 ENSMUSG00000006205 202. -0.294 0.198 -1.43 0.153 1.00
## 4 ENSMUSG00000020212 33.8 0.134 0.192 0.633 0.527 1.00
## 5 ENSMUSG00000020541 745. 0.0694 0.148 0.449 0.654 1.00
## 6 ENSMUSG00000020719 6565. -0.00388 0.126 -0.0323 0.974 1.00
## 7 ENSMUSG00000020810 125. -0.124 0.193 -0.534 0.593 1.00
## 8 ENSMUSG00000020865 681. -0.0432 0.186 -0.230 0.818 1.00
## 9 ENSMUSG00000021356 99.6 -0.133 0.193 -0.744 0.457 1.00
## 10 ENSMUSG00000021428 142. -0.0579 0.156 -0.372 0.710 1.00
## # … with 39 more rowspc2_dif## # A tibble: 49 x 7
## X1 baseMean log2FoldChange lfcSE stat pvalue padj
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ENSMUSG00000000441 1066. -0.0407 0.155 -0.269 0.788 1.00
## 2 ENSMUSG00000001870 179. 0.0645 0.198 0.552 0.581 1.00
## 3 ENSMUSG00000008575 523. -0.158 0.196 -0.725 0.468 1.00
## 4 ENSMUSG00000009555 612. 0.216 0.165 1.29 0.198 1.00
## 5 ENSMUSG00000010392 679. 0.00609 0.144 0.0390 0.969 1.00
## 6 ENSMUSG00000012350 1293. 0.106 0.195 0.579 0.563 1.00
## 7 ENSMUSG00000017418 377. -0.149 0.183 -0.827 0.408 1.00
## 8 ENSMUSG00000017478 698. 0.0420 0.196 0.227 0.821 1.00
## 9 ENSMUSG00000019960 736. 0.0174 0.185 0.0831 0.934 1.00
## 10 ENSMUSG00000020149 4293. -0.0580 0.112 -0.531 0.596 1.00
## # … with 39 more rowspc1_DEG <- subset(pc1_dif, pc1_dif$padj < 0.1 & abs(pc1_dif$log2FoldChange > 0.1))
pc2_DEG <- subset(pc2_dif, pc2_dif$padj < 0.1 & abs(pc2_dif$log2FoldChange > 0.1))
pc1_genes_DEG <- pc1_genes[pc1_genes$X1 %in% pc1_DEG$X1,]
pc2_genes_DEG <- pc2_genes[pc2_genes$X1 %in% pc2_DEG$X1,]
pc1_genes_DEG <- inner_join(pc1_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
pc2_genes_DEG <- inner_join(pc2_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
print(paste("The number of MK2 inhibitor targets from TNF model from UC PC1:", dim(pc1_genes_DEG)[1]))## [1] "The number of MK2 inhibitor targets from TNF model from UC PC1: 0"print(paste("The number of MK2 inhibitor targets from TNF model from UC PC2:", dim(pc2_genes_DEG)[1]))## [1] "The number of MK2 inhibitor targets from TNF model from UC PC2: 0"